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ABSTRACT, 


Recently,  Saunders  1974  has  general  i zeathe  sir- cal  led  reciprocal  property 
for  normality.  Distributions  having  this  property  are  called  ^-normal , and 
it  is  of  interest  to  make  statistical  inferences  about  the  relevant  parameters 
when  sampling  from  such  a distribution.  Previous  work  on  such  problems  has 
been  from  the  sampling  viewpoint. 


In  this  paper,  we  approach  the  inference  problem  from  the  Bayesian  point 

of  view  and  investigate  the  posterior  of  the  parameters  involved  when  sampling 

is  from  the  c-normal  with  parameters  a and  8.  Two  special  cases  are 

examined,  £(v)  = loggv,  which  gives  rise  to  the  lognormal  distribution,  and 
1 /2  -1 /2 

S(v)  = v ' - v , a case  that  arises  in  certain  fatigue  problems  (Saunders 

and  Birnbaum  lQfiQ) , — — — ■ 
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SIGNIFICANCE  AND  EXPLANATION 


' 


Statistical  problems  involving  "time  to  failure",  "strength  of  material", 
and  "fatigue"  are  often  characterized  by  a statistical  distribution  known  as 
the  lognormal  distribution;  the  length  of  life  of  a system  subjected  to  fluctuat- 
ing stresses  by  a periodic  loading  can  sometimes  be  predicted  through  the  use 
of  the  Birnbaum-Saunders  life  length  distribution. 

These  distributions  are  special  cases  of  a more  general  class  of  distribu- 
tions based  on  certain  functions  s(t),  that  map  0 £ t £ °°  into 
-«  £ c(t)  £ <=°  and  have  the  property  that  s(t)  = -c(l/t).  In  the  above  cases, 
we  have  s(t)  = an  t and  s(t)  = t1^2  - t-1/^2,  respectively. 
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BAYESIAN  INFERENCE  ON  PARAMETERS  ASSOCIATED  WITH  FAMILIES 

* 


CLOSED  UNDER  RECIPROCATION 
Norman  R.  Draper  and  Irwin  Guttman 


1.  INTRODUCTION 


** 


Recent  studies  on  reliability  analysis  in  fatigue  testing  have  focused 
on  random  variables  distributed  under  the  5-normal  distribution;  see,  for  ex- 
ample, Birnbaum  and  Saunders  1968,  1969;  Esary  and  Marshall  1970;  Marshall  and 
Proschan  1972;  Saunders  and  Birnbaum  1969;  and  Whittaker  and  Besuner  1969.  A 
precise  definition,  given  by  Saunders  1974  may  be  stated  as  follows. 

Definition  1.  Suppose  5 is  a concave, monotone  increasing*  differentiable  function 
from  (C,»)  onto  (-»,°°),  and  such  that 


£(t)  = -C(l/t),  t > 0. 


(1.1) 


Then  the  non-negative  random  variable  T is  said  to  have  the  5-normal  dis- 
tribution with  parameters  (a,B),  where  a > 0,  6 > 0,  if 


jC(T/3)  = Z, 


(1.2) 


where  Z Is  a standard  normal  random  variable,  that  is  Z - N(0,1). 
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An  equivalent  definition  statement  is  that  the  non-negative  random 
variable  T is  5-normal  if 

T = 3<|»(a7),  a. 3 > 0,  (1.3) 

where  <>  * 5 ^ is  a convex  monotone  increasing  transformation  from  (-«,»)  onto 
(0,*»)  with  the  property  that 


*(-v)  = l/t|»(v)  (1.4) 

and  where  Z - N(0,1).  It  is  convenient  to  use  (1.3)  to  establish  the  first 
part  of  the  following  theorem,  due  to  Birnbaum  and  Saunders  1969  and  Saunders  1974. 
The  second  part  is  obvious. 

Theorem  1.  If  T is  5-normal  with  parameters  a,B  > 0,  then  1/T  is  5-normal  with 
parameters  a,  1/3.  Furthermore,  for  any  real  a > 0,  aT  is  5-normal  with  parameters 
a,a0. 

To  quote  Saunders  1974,  p.  534,  the  above  shows  that  "...  the  parameter  3 
Is,  in  fact,  a scale  factor  for  any  5-normal  distribution."  Usinj  (1.1),  it  is 
quickly  verified  that  3 is  "...  also  the  median.  It  is  called  the  characteristic 
life  in  many  applications".  To  see  that  3 is  the  median,  we  note  from  (1.3) 
and  the  fact  that  * 5~\  that 

P(T  < 3)  - P(Z  < cf^O));  (1.5) 

However,  it  Is  obvious  from  (1.1)  that  5(1)  ■ 0,  whereupon  the  result  follows. 

In  this  paper,  we  discuss  the  Bayesian  approach  to  inferences  about  a and 
0.  In  Section  2,  we  obtain  the  joint  posterior  distribution  of  (a, 3),  and  the 
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The  former  case  gives  rise  to  the  lognormal  distribution,  one  often  assumed  in 
reliability  work.  The  latter  case  stems  from  a model  for  fatigue  failure  given 

by  Saunders  and  Birnbaum  1969.  It  is  easily  verified  that  (i)  and  (ii)  satisfy  (1.1) 
and  that  their  inverses  satisfy  (1.4).  Finally,  in  Section  4,  we  discuss  the  question 
of  a Bayesian  tolerance  interval  for  B,  which  has  applications  to  the  question  of 
warranty  life  of  a system  involving  k loading  spectra. 


2.  THE  POSTERIOR  OF  (a, a) 


The  probability  density  function  (pdf)  of  Z,  the  standard  normal 
variable,  is 


♦(2)  = (2ir)~1/2  exp  {-z2/2}.  (2.1) 

The  Jacobian  of  the  transformation  (1.3)  is,  from  (1.2), 

jf"  (C'(t/6))/(a6),  (2.2) 

where  the  prime  denotes  differentiation  with  respect  to  t.  This  derivative  is 
positive  since  K is  increasing.  Hence  the  pdf  of  T,  given  a, 6,  is 

p(t|a,B)  * *[a"?£(t/B)]S‘(t/e)/(ae)  (2.3) 

if  t > 0,  and  zero  otherwise.  Note  that,  in  view  of  (2.1),  (2.3),  and  (1.3), 
we  may  express  the  mean  and  variance  of  the  distribution  (2.3)  as 


E(t|a,3)  - / t ♦[o'15(t/B)]C(t/B)/(aB)dt 
0 

(2.4a) 


■ / 3</»(az)  $(z)dz, 


and 
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(2.4b) 


Var(t|a,B)  = / t2  ♦[a‘1C(t/B)]C (t/B)/(aB)dt  - [E(t|a,8)]2 
0 


- / BV(az)<j>(z)d2  - [E(t|a,B)]2, 


respectively.  That  is,  (2.4a)  and  (2.4b)  are  functionals  of  £ (or,  equivalently 
il>),  and  can  be  evaluated  for  specific  5 (or  4»);  see  Section  3 for  examples. 

Suppose  now  that,  prior  to  taking  observations  on  the  random  variable  T, 
whose  pdf  is  given  by  (2.3),  very  little  is  known  about  the  parameters  a and  B. 
In  this  situation,  we  can  use  the  Jeffreys  invariant  prior 

p(o,B)  « 1 1 (cx,e)  | 1/2  (2.5) 

where  I Is  the  usual  information  matrix  consisting  of  negative  expected  (with 
respect  to  t)  second  derivatives  (with  respect  to  a and  or  B)  of  Anp(t|a,B) . 

The  likelihood  function  l of  (a, 8)  for  n observations  is 


n n 1 

£(a,B|t) * (aB)  " ir  (♦[a',C(t1/B)]C'(t1/B)}.  (2.6) 

1=1  1 1 

Combining  (2.5)  with  (2.6)  via  Bayes'  Theorem  leads  to  the  posterior  of  (a, 6), 
given  by 

p(a,B|t)  « (aB)"n[  it  ♦[a‘1C(ti/B)]5,(ti/B)]  |I(a,B)|1/2  (2.7) 

1=1  1 

for  a,B  > 0,  and  zero  otherwise.  In  general,  the  marginal  posterior  distributions 
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of  a and  of  8 are  obtained  by  integrating  out  the  unwanted  parameter. 
In  cases  where  a and  3 are  independp"*,  a priori , so  that 


|I(a,8)|1/2  « h1(o)h2(B).  (2.8) 

a situation  that  occurs  in  the  two  examples  we  provide  later, 
we  may  express  the  marginal  (posterior)  distribution  of  the  median  8,  the 
characteristic  life,  as 


n od 

p(e|t)  « B’V(B)  [ TT  5'(t./e)]  / a~nh, (a)  exp{-52(t./g)/2a2}da,  (2.9) 

^ 1=1  1 0 1 1 

and  simllarily,  the  marginal  (posterior)  of  a becomes 


oo  n 

p(a|t)  « a"V(a)  / e"V(g)  [ * {♦[a"15(tl/B)]€'  (t,/8)}d8.  (2.10) 

1 0 ^ 1=1  1 1 

It  is  interesting  to  note  from  (2.7)  that  the  parameter  8 has  a marginal 

(posterior)  distribution  that  exhibits  a reciprocal  property.  For  if  we  set 

-1  -2 
8 = 6,  with  Jacobian  38/30  = -0  , and  use  the  facts  that 


02  |I(a.9)|1/2  = |I(a,8)|1/2, 
<fr(z)  = <p(-z), 
c(t/e)  - -e(e/t), 
t2c* (t/8)  » 82C'(B/t), 


(2.11) 
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we  obtain  a joint  distribution  p(a,e|r)  with  r = (r^rg,...,!*  )'  and  ri  = l/ti 
of  exactly  the  same  form  as  (2.7).  In  other  words  the  posterior  distributions 
of  0 and  of  0 ^ belong  to  the  same  family,  a property  that  arises  from  the 
assumption  of  5-normality. 

Tolerance  regions  for  future  life-lengths. 

We  here  review  standard  definitions;  see  Guttman  (1967). 

Definition  2:  Suppose  (tj ,t2>. . . .tn)  = t'  is  a vector  of  n independent  ob- 
servations sampled  from  the  distribution  p( t [ 0 ) . Then  a region  of  t values 
defined  by  a function  S(t-| ,t2«. . . ,tn) CR1  = {t|-°°  < t < »}  is  said  to  have 
‘ coverage  C if 

C = C(S)  = / p(t|0)dt.  (2.12) 

S 

Note  that  C,  a functional  of  S,  is  a function  of  6. 

Definition  3:  S is  a 6-expectation  tolerance  region  if 

E{C(S) |t)  * 6 (2.13) 

the  expectation  being  taken  with  respect  to  the  posterior  distribution  p(e|t). 

A procedure  for  the  construction  of  an  S satisfying  (2.13)  is  as  follows. 
(For  proofs  and  discussion,  see  Guttman  1967;  1970,  pp.  132-3.) 

1.  Find  the  predictive  distribution 

h(t|t)  = I p(t|e)p(e|t)de  = / p(e,t| t)de  (2.14) 

ecft  ~ ~ ~ een  ~ ~ ~ 

where  ft  denotes  the  parameter  space. 
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2.  Find  S,  a 6-confidence  region  for  t such  that 


P(tr^|t)  « / h(t|t)dt  = 6 . (2.15) 

S 

For  our  problem,  with  9 = (a,8)\  p(t|0)  Is  (2.3)  and  p(0|t)  Is  (2.7)  assuming 

the  prior  of  (2.5).  It  follows  that  p(0,t|t)  is  of  the  form  of  (2.7)  with 

(n+1)  replacing  n and  the  vector  t extended  to  (t^.t,,...,!  ,t)  * t'  ,, 

- i c n _n+i 

say.  Equation  (2.14)  then  defines  the  normalizing  constant  of  this  distribution. 
This  process  will  be  Illustrated  in  Section  3. 


3.  TWO  SPECIAL  CASES 


r "”f 


The  Log-Normal  Distribution. 

In  many  reliability  papers,  the  distribution  of  the  random  variable  being 
studied  (e.g.,  time  to  failure,  strength,  fatigue)  is  assumed  to  be  log-normal. 
That  is, 

£(v)  = Jlnv  (3.1) 

so  that,  in  the  notation  of  Definition  1, 

a_1£n(T/e)  = Z,  or  T = 6 exp{aZ},  (3.2) 

where  Z - N(0,1).  If  we  apply  the  foregoing  work,  we  find  a number  of  well- 

-2  -1 

known  results.  The  invariant  prior  is  proportional  to  a 8 and,  a posteriori, 
W/ot  ~ xn  where  W = E(u^-u)  and  u^  = tnt^,  and  n(f,n0-u)/W  - tn-  We  omit  the 
details;  for  closely  related  results  with  a different  prior,  and  for  a number 
of  other  developments,  see  Zellner  (1971).  For  comments  on  the  improper  prior 
distribution,  see  the  second  example. 


The  Birnbaum-Saunders  Life  Length  Distribution. 

Birnbaum  and  Saunders  1969  are  concerned  with  a random  variable  T,  the  length 
of  life  of  a system  which  is  subjected  to  fluctuating  stresses  by  a periodic  loading. 
Under  various  assumptions,  they  show  that,  in  the  language  of  Definition  1,  T is 
C-normal  with 


vV2  . v-1/2| 


(3.3) 


that  is,  T is  such  that 


o'Ud/B)  = cf'cd/e)172  - (fJ/T)"1/2]  = Z. 


It  is  straightforward,  but  tedious,  to  verify  that,  for  (3.4) 


|e[(aZ)2+2  - {(aZ)4+4(aZ)2}1/2]  if  Z < 0, 
]fi[(aZ)2+2  + {(aZ)4+4(aZ)2}1/2  if  Z > 0, 


so  that  (see  (1.3)  with  W * aZ), 


i(W2+2  - {W4+4W2}172]  if  W < 0, 
J{W2+2  ♦ (W4+4W2)1/2j  if  W > 0. 


(3.4) 


(3.5) 


(3.6) 
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Is  convex,  monotone  increasing,  and  such  that 


*(W)  = l/i|»(-W). 

Using  properties  of  Z = N(0,1)  we  find,  from  (3.5),  that 
E(T|o,B)  = $ {l+]a2}, 

E(T2|a,0)  = 02{l+2a2+3a4/2), 


(3.7) 


(3.8) 


so  that 


Var(T|a,0)  = 02a2{l+5a2/4), 


(3.9) 


with 

{Var(T|a,0)}1/2/E(T|a,0)  = a(l+5a2/4)1/2/(l+a2/2),  (3.10) 

• function  of  a alone.  The  results  (3.8),  (3.9)  were  found  by  Birnbaum  and 
Saunders  1969,  using  a different  method.  Of  course.  If  we  substitute  (3.4)  in 
(2.3),  we  find  inmediately  that 

p(t|a,0)  - {(t+8)/[a207r(2t)3]1/2}  exp{-[(t/B)1/2-(8/t)1/2]2/(2a2)}  (3.11) 


for  t > 0,  because 


5'(t/B)  = 61/2(t+8)/(2t3/2} 


(3.12) 


and  we  may  find  moments  directly  from  (3.5). 

The  Jeffrey's  invariant  prior  for  this  problem  has  the  interesting  form 
(details  are  given  in  the  Appendix): 


p(o.B)  - |I(a,8)|1/2  = [a2(l-G(a))+2]1/2/(a2B) 


(3.13) 


where 


G(a)  = (l-4>(2/a)}/ (a<j)(2/a)}. 


with  4>(u)  = f 4>(z)dz.  If  we  set  x = 2/a,  (3.14)  takes  the  form 


(3.14) 


G(2/x)  = |x{l-4>(x)}/4>(x)  = ^xR(x) 


(3.15) 


where  R(x)  is  Hill's  ratio  as  given  by  Johnson  and  Kotz  1970,  p.  278,  (72). 

The  prior  (3.13)  is  improper,  that  is,  the  integral  with  respect  to  a 
and  8 diverges;  this  Is  true  with  respect  to  the  parameters  individually,  in 
fact.  In  most  practical  cases,  this  is  not  a matter  for  concern.  Improper  priors 


often  arise  in  the  expression  of  imprecise  prior  knowledge,  and  correspond  to 

"flat"  prior  information  in  some  suitably  transformed  parameter  space.  Once 

the  information  in  the  likelihood  is  factored  in,  the  consequent  posterior 

typically  converges.  If  it  did  not,  the  sample  information  would  be  extremely 

weak  and  of  the  same  character  as  the  prior  information.  Improper  priors 
can  be  regarded  as  approximations  to  proper  priors  in  situations  where  the 

proper  prior  is  "flat"  in  the  region  covered  by  the  likelihood.  In  situations 

where  a Bayesian  solution  is  paralleled  by  the  sampling  theory  solution,  improper 

priors  frequently  occur  in  this  way.  In  such  cases,  criticism  of  an  improper 

prior  solution  implies  criticism  of  the  sampling  theory  solution  as  well. 

For  recent  discussion  on  the  various  aspects  of  improper  priors,  see  Dawid, 

Stone  and  Zidek  (1973). 

Suppose  now  that  n independent  observations  (t-j  ,tg tn)  = t'  are 

taken  on  T,  whose  distribution  is  (3l 1 1 ) . From  (2.7),  (3.12)  and  (3.13),  and 
after  some  simple  algebra,  the  posterior  of  (a,B),  given  t,  is  seen  to  be 


p(a,B|t)  = Ka'(n+2)B'(n+2)/2{a2[l-G(a)>2}1/2[  " (6+t-)] 

i=l  1 

X exp{-[S2(B-e0)2+V]/[2a2B]}, 


(3.16) 


where  K is  a normalizing  constant. 


I 


(3.17) 

(3.18) 

Using  the  Schwarz  inequality,  it  is  easily  shown  that,  because  t^  > 0,  V ^ 0. 

For  this  case,  the  marginal  distributions  p(a|t)  and  p(8|t)  are  not 


V = Srn2/S2  = (S1S2-n2)/S2, 


n n _i 

S,  = E t. , S9  = E tT  , and  = n/S9 
1 1=1  1 z 1=1  1 0 2 
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convenient  explicit  forms  and  so,  in  working  our  illustrative  example  below, 
we  have  integrated  (3.16)  numerically  over  6 and  a respectively.  (We 
could,  of  course,  express  p(a|t)  in  series  form  explicitly.) 

Numerical  Illustration 

We  used  the  Rand  Corporation  Tables  1955  to  generate  15  random  z^'s  which 
were  converted  into  t^'s  via  (3.5).  These  are  given  in  Table  3.1;  it  follows 
that  n * 15,  V = 454.39,  Sj  s 498.50  Sg  * 5.10,  and  0Q  = 2.94,  to  two  decimal  places. 

Table  1.  Fifteen  independent  observations  of  life  length. 

45.3  35.8  115.7  99.9  52.5  66.9  42.4  13.3 

7.1  0.4  1.8  8.5  1.0  2.5  5.4 

The  posterior  distributions  for  a and  6,  obtained  by  numerical  integration  of 
(3.16)  with  respect  to  8 and  a,  respectively,  are  shown  in  Figures  1(a)  and  (b). 
From  this  calculation,  we  obtain  95%  posterior  confidence  intervals  for  a and 

6 as 


1.73  < a < 3.38  4.57  < 6 < 16.59 


the  true  values  being  2 and  10  respectively. 

The  mean  and  mode  of  p(a|t)  are  at  2.41  and  about  2.17,  while  the  mean  and 
mode  of  p(B|t)  are  at  9.53  and  about  8.02,  respectively. 

The  question  of  suitable  approximations  for  the  marginals  remains  to  be 
answered. 


For  this  example,  h(t|t)  of  (2.14)  is,  by  extension  of  (3.43),  proportional 
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Posterior  distributions  of  o and  3 for  the  life  length  numerical 


0 90.0 


to  the  reciprocal  of 


/ / a"^n+3^8‘^n+3^/2{a2[l-G(a)]+2}1/2/(t+e)t_3/2  ! 
0 0 1=1 

exp{-[T2(8-3;)2+Tr(n+l)2/T2]/(2a28)}dad8, 


(ti+B)tT3/2  x 


(3.19) 


where  T,  = t+(t,+t2+. . .+tn)  = t+Sr  T2  = t‘1+(t"1+t21+. . ,+tjJ1 ) = t"1  + S? 
and  8^  = (n+l)/T2  = 8Qt(n+l)/(n+80).  There  is  no  closed  form  for  (3.19)  and 
so,  for  a given  set  of  data,  it  must  be  evaluated  numerically  for  a selected 
grid  of  t values  and  then  normalized  by  a third  numerical  integration.  During 
the  latter  integration  the  various  percentage  points  (e.g.  0.025  and  0.975) 
which  are  the  ends  of  the  corresponding  tolerance  regions  (e.g.  95%)  can  also  be 
evaluated  and  recorded. 

We  did  this  for  the  data  in  Table  1 using  the  supplementary  figures 
recorded  just  above  that  table,  and  obtained  the  following  tolerance  regions  for 
a future  observation:  (1.71,87.50)  for  a 0.90-expectation  region,  (1.10,  93.47) 
for  a 0.95-expectation  region  and  (0.52,  98.64)  for  a 0.99-expectation  region. 
Comparison  with  Saunders*  estimates. 

Saunders  (1974)  offers  a variety  of  alternative  estimates  for  8 and  a. 


from  the  third  display  on  p.  535 


a*  = 5.855,  8*  = 11.343; 

from  p.  535  in  Remark  3.6, 


B - 12.280; 


and  from  p.  536,  second  display 


B°  = 13,3,  the  middle  Tj. 

A A 

We  see  that  only  a,  8.  and  8 provide  values  that  appear  appealing  compared 
with  those  we  obtain  from  our  Bayesian  approach,  and  from  our  knowledge  of 
the  true  values  o = 2,  8 * 10. 
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APPENDIX 


Jeffrey's  invariant  prior  p(g,B)  for  Section  3.  second  example. 
From  (3.11)  we  have. 


L = £np  = constant  + £n(t+$)  - £na  - ^£nB  - (t8"^+t~^8-2)/(2a2).  (A.l) 


The  expected  (with  respect  to  t)  matrix  of  negative  second  derivatives  is 


32L 

32LH 

= E 

3a‘4(tB-1+t'1B-2) 

-a-2  a'3(tB'2-t'1) 

3a^ 

3a3B 

32L 

32L 

- o 

a”3(tB"2-t~1 ) 

(t+B)"2+ta'2B"3-Je'2 

3a3B 

3^ 

L 

c 

s 

!11  *12 

9 

( 

_*12  *22_ 

(A. 2) 


say.  Now 


E(3(£np)/3a)  = 0 implies  E(tB”^+t_1B-2)  = a2. 


(A. 3) 


It  immediately  follows  that  = 2/a  . Also 


T 


i 


1 1 
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(A. 


E(3(£np)/3B)  = 0 implies  E(t+0)-1  - ^E(t  -t0~Z)/aZ  = |e-1 

and,  after  some  algebra,  it  can  also  be  shown  that 

EU+fS)'1  = ^B’1.  (A.! 

It  follows  from  (A. 2),  (A. 4)  and  (A. 5)  that  = 0.  Now  we  can  show  directly 
that 


E(t)  = |e(2  + aZ), 


(A.! 


whereupon,  from  (A. 2), 


I22  = E(t+B)'Z  + a"VZ.  (A.( 

-2 

We  now  need  E(t+0)  . This  is  a tedious  calculation  in  which  account  must  be 

taken  of  (3.5),  and  a series  of  reductions  performed.  Eventually  we  find  that 

E(t+0)"2  = 1b'2{1-G(cx)}  (A.5 

where  G(a)  is  defined  by  (3.14).  Combining  (A. 2),  (A. 6)  and  (A. 7)  and  the 
facts  that  ■ 2a‘Z,  I.?  = 0 now  provides  the  result  (3.13). 
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